%% On the Optimal Design of Transfers and Income-Tax Progressivity
% Read results from "info" file for optimization grid searches

clear; clc; close all;

%% Read results

% Import "info"
info = readmatrix('info.txt');      % read info matrix
info = info(any(info,2),:);         % remove zero line
ninfo = size(info,1);               % number of tax systems

% Replace transition welfare if non-convergence
for i = 1:ninfo
    if (info(i,8) == 0 || info(i,5) > 1e-5 || info(i,6) > 1e-5)
        info(i,8) = 10000;
    end
end

% Sort
info_sorted_trans = sortrows(info,8);           % by welfare including transitions
info_sorted_m_theta = sortrows(info,[3,2]);     % by m and theta

% grid for m
mgrid = unique(info_sorted_m_theta(:,3));       % grid for m 
nm = size(mgrid,1);                             % number of points

% values for theta
theta_values = unique(info_sorted_m_theta(:,2));% unique theta values
ntheta_max = size(theta_values,1);              % number of points
theta_welf = zeros(ntheta_max,6);               % auxiliary theta and welfare storage
theta_grid = 0*mgrid;                           % grid to store optimal theta 
transfer_grid = 0*mgrid;
amtr_diff_grid = 0*mgrid;
mtr_top10_grid = 0*mgrid;
amtr_grid = 0*mgrid;

% Find optimal thetas for each m
for i = 1:nm
    m = mgrid(i);                                               % pick m
    iter = 1;                                                   % auxiliary counter                                                  
    theta_welf(:,2) = 10000;                                    % set 
    for j = 1:ninfo
        if (info_sorted_m_theta(j,3) == m)                      % consider only elements of "info" with the right m
            theta_welf(iter,1) = info_sorted_m_theta(j,2);      % pick theta
            theta_welf(iter,2) = info_sorted_m_theta(j,8);      % pick corresponding welfare
            theta_welf(iter,3) = info_sorted_m_theta(j,17);
            theta_welf(iter,4) = info_sorted_m_theta(j,19)-info_sorted_m_theta(j,20);
            theta_welf(iter,5) = info_sorted_m_theta(j,19);
            theta_welf(iter,6) = info_sorted_m_theta(j,20);
            iter = iter+1;                                      % add to auxiliary counter
        end
    end
    [~,loc] = min(theta_welf(:,2));                             % pick best of all the thetas computed for this m
    theta_grid(i) = theta_welf(loc,1);                          % value of best theta for this m
    transfer_grid(i) = theta_welf(loc,3);
    amtr_diff_grid(i) = theta_welf(loc,4);
    mtr_top10_grid(i) = theta_welf(loc,5);
    amtr_grid(i) = theta_welf(loc,6);
end

% Find best m-theta combination
[~,opt_loc] = min(info_sorted_m_theta(:,8));
theta_opt = info_sorted_m_theta(opt_loc,2);
m_opt = info_sorted_m_theta(opt_loc,3);

% Save mgrid and rescaled optimal thetas
m_grid_bench = mgrid;
theta_grid_bench = theta_grid/2;
transfer_grid_bench = transfer_grid;
amtr_diff_grid_bench = amtr_diff_grid;

% plot(transfer_grid/1000,amtr_diff_grid,transfer_grid/1000,mtr_top10_grid,transfer_grid/1000,amtr_grid,'Linewidth',2)
% legend('AMTR Top 10 - AMTR All','AMTR Top 10','AMTR All','Location','best')
% xlabel('Transfer in Thousands of Dollars')

filename='mtheta_bench.mat';
save([filename],'transfer_grid_bench','amtr_diff_grid_bench','m_grid_bench','theta_grid_bench') 
